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Abstract 

A numerical method is developed leading to Lyapunov operators to approximate the 
solution of two-dimensional Boussinesq equation. It consists of an order reduction 
method and a finite difference discretization. It is proved to be uniquely solvable and 
analyzed for local truncation error for consistency. The stability is checked by using 

Lyapunov criterion and the convergence is Some numerical implementations 

are provided at the end of the paper to validate the theoretical results. 
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1 Introduction 

The present paper is devoted to the development of a numerical method based 
on two-dimensional finite difference scheme to approximate the solution of the 
nonlinear Boussinesq equation in M? written on the form 

u tt = Au + u xxxx + (u 2 ) xx , ((x,y),t) e O x (t ,+oo) (1) 

with initial conditions 

u 

u(x, y, t ) = uo(x, y) and -wr(x, y, t ) = <p(x, y), (x, y) G Q (2) 

at 
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and boundary conditions 



du 
dn 



(x,y,t) = 0, ((x,y),t) G dn x (to, +00) 



(3) 



on a rectangular domain Q =]Lq, Li[x]Lo, L%[ in R 2 . to > is a real parameter 

fixed as the initial time, u u is the second order partial derivative in time, 

A = ^2 + -§^2 is the Laplace operator in IR 2 , u xx and m iixi are respectively 

V d 
the second order and the fourth order partial derivative according to x. — — 

on 

is the outward normal derivative operator along the boundary dVl. Finally, u, 
uq and ip are real valued functions with uq and (p are C 2 on Q. 

The Boussinesq equation 



Our main idea consists in developing a numerical scheme to approximate 
the solution of ([I])-®. The time and the space partial derivatives are re- 
placed by finite-difference approximations in order to transform the initial 
boundary- value problem ([I])-© into a linear algebraic system. An order re- 
duction method is adapted leading to a system of coupled PDEs which is 
transformed by the next to a discrete algebraic system. The resulting method 
is analyzed for local truncation error and stability and the scheme is proved 
to be uniquely solvable and convergent. 



2 Two-Dimensional Finite Difference Scheme 

Consider the rectangular domain Q =]Lq,Li[x]Lq,Li[g M 2 and an integer 
J G N*. Denote h = — - for the space step, Xj = L + jh and y m = 

rj 

Lq + mh for all (j, m) e I 2 = {0, 1, . . . , J} 2 . Let I = At be the time step and 
t n = to + nl, n G N for the discrete time grid. For (j, m) G I and n > 0, u™ m k 
will be the net function u(xj,y m ,t n ) and U? m the numerical solution. The 
following discrete approximations will be applied for the different differential 
operators involved in the problem. For time derivatives, we set 



jjn+l _ jjr, 
3,™ ]■ 



rn— 1 



and u u ~^ 



m+l 



- 2U n 




+ m 



jn—1 



21 



and for space derivatives, we shall use 



for first order derivatives and 

U XX ~H — , Uyy — 

for second order ones, where for n£ff* and a£l, 

u n ' a = a U n+1 + (1 - 2a)U n + a IP' 1 . 

P 1 
Finally, we denote a = — and 8 — —. 

h 1 h 1 



3 Discrete Two-Dimensional Boussinesq Equation 



The method is based on an order reduction idea in which we set 

v = u xx + u 2 . (4) 

For (j, m) G I 2 an interior point of the grid J 2 , (J = {1, 2, . . . , J — 1}), and 
n > 1, the following discrete equation is deduced from (pQ). 

Tjn+l _ 2Jjn , Tjn-1 

— era ( JJ n+1 — 4/7 n+1 4- JJ n+1 4- JJ n+l 4- JJ n+1 \ 

+a(l - 2a) ([/;_ lim - 4E/£ m + l7J» +lffB + E/*^ + f/« m+1 ) 

+aa (tfj*~£ m - 4[/™~ 1 + Uj+l m + Uf^_ x + U?m+i) ( 5 ) 

-i-rrni ( V n+1 — 9V n+1 -4- 
T-fJO! \^'j—ljrn ZV j,m ' V j+l,n 



+a(l - 2a) (VJl l>m - 2V£ m + V? +lt „ 

+aa (y^ m - 2V%? + v%r£ m ) • 

Similarly, the following discrete equation is obtained from equation 
V^ 1 + V?- 1 = 25a (Vp* n - 2XT£t + U?+i,m) 

rn — OTjn _i_ Tin 
j—l,m ^j,m 1 w j'+l,m 



+26(l-2a)[U?_ hm -2UZ m + U. 
+25a (V 



+2Sa (U£l n - 2U?- 1 + U£l n 



(6) 



pn-l _|_ pn 

where F(u) = u 2 , F n = F(u n ) and F n = . The discrete boundary 

conditions are written for n > as 

U? m = U« and U«_ l>rn = U« (7) 



3 



E£i = E/£_i and U^_, 



U 



Next, as it is motioned in the introduction, the idea consists in applying Lya- 
punov operators to approximate the solution of the continuous problem ([I])-© 
or its discrete equivalent system flS])-®. Denote 



ai = — + 2aa, a 2 



-acr, 



bi = 1 - 2(1 - 2a)a, b 2 = (1 - 2a)a, 
c\ — (1 — 2a)<5 and c 2 = ct<5. 



Equation ([5]) becomes 



+a 2 (yj!^ m — 2V™+ 1 + V^i| TO J 

^Uf_i^ m + &iC^ m + b 2 Uf +lm + b 2 U™ m _ 1 + biU™ m + b 2 U? m+1 



—a 2 Uj_ lm 



-a 2 ( 



n— 1 



j — l,m 



2K 



n— 1 



n— 1 



(9) 



Equation ([6]) becomes 



v j,m zc 2 ^j-l,m ZU i,m ' u j+l,m 

2ci {u r j_ l m — 2U™ m + ^JYi,™) 
+2c 2 (Uj-i tm — 2f/ ? " m 1 + t/™ + i;m) 
-W^ 1 + 2F(tHf 



Denote A, I? and it! the matrices defined by 

An 2a 2 \ 

d 2 <X\ d 2 ' ■ ' ■ '. 

o '•. '•. '-. '•. ; 

: ■•. "•• ■•• ■•. 

'. ' • ' • a 2 (i\ Qj 2 

\0 2a 2 at J 



A 



B 



(b x 2b 2 \ 

b 2 bi b 2 : 

o '-. '•. '•. '•. ; 

: ■•• "•• "•• ••• 

: b 2 bi b 2 

\0 26 2 bij 



(10) 
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R 



and 

(-2 2 \ 

1 -2 1 

o '•. '•. '•. '•. ; 

; ■•. ■•. ■•. ■•. o 

: 1 -2 1 

\0 2 -2 J 

The system (17|)- (TTU|) can be written on the matrix form 

C A {U n+1 ) + a 2 RV n+1 = C B {U n ) - Ca^- 1 ) + R(b 2 V n - a 2 V n ~ l ) 

-V n+1 - c 2 RU n+1 = R( Cl U n + c 2 XJ n - x ) - -V^ + F^ 
2 2 



'id 



for all n > 1 where 

U n = (U2 m )0<j,m< J-i V 1 — {V£m)0<3,m<J an d = (F(U£ m ))o<j,m<J 

and for a matrix Q e Cq is the Lyapunov operator defined by 

C Q (X) = QX + XQ, VI e M {J+ir (R). 



4 Solvability of the Discrete Problem 

We are now ready to set our first main result. 

Theorem 1 The system ( TO]) is uniquely solvable whenever U° and U 1 are 
known. 

Proof. It reposes on the inverse of Lyapunov operators. Consider the endomor- 
phism $ which associates to any element (X, Y) in .M(j+i)2(R) x AArj+xy(M) 

its image Y) = (AX + XA + a 2 RY, W - c 2 RX). To prove Theorem QJ 
it suffices to show that /cer$ is reduced to 0. Indeed, 

Y) = {AX + XA + a 2 RY, l -Y - c 2 RX) = (0, 0) 

or equivalently, 

Y = 2c 2 RX and (A + 2a 2 c 2 R 2 )X + XA = 0. 
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So, the problem is transformed to the resolution of a Lyapunov type equation 
of the form 

C WA (X) = WX + XA = (12) 
where W is the matrix given by W = A + 2a2C2-R 2 . Denoting 

co = 2a 2 C2, co\ — ai + 6co, coi = co\ + co and co 2 = a<i — 4o> 

the matrix is explicitly given by 

^ co\ 2co 2 2co 

0J2 oJ\ co 2 co ' • • ' • • ' • • : 

co U2 uji 0J2 co : 

'•• '•• '•• '•• '•• '•• : 



W 







V 



Next, we use the following result. 



CO CO2 COi CO2 CO 
' ■ ■ CO CO2 CO 1 CO2 

. . . 2co I002 001 ^ 



Lemma 4.1 Let E be a finite dimensional (R orC) vector space and ($ n ) n be 
a sequence of endomorphisms converging uniformly to an invertible endomor- 
phism $. Then, there exists n such that, for any n > n , the endomorphism 
$ n is invertible. 

Assume now that / = o(h 2+s ), with s > which is always possible. Then, the 
coefficients appearing in A and W will satisfy as h — > the following. 



-+eh 2+2s 



For 1 < % < J -I, 

A^i-i 

For 2 < i < J -2, 
Similarly, 



Ai 



i+i 



2 2 



0. 



W ,o = Wj,j = - + 2ae h 2+2s - I2a 2 e ti 



2s 



1 

2" 



W ltl = Wj^j^ = - + 2ae h 2+2s - Ua 2 e h 2s 



1 

2 
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and 



Li— 1 



Wi 



W A Wjr, 



Finally, 



W, 



t,t-2 



2 2 



aeh 2+2s + 8a 2 sh 2s 



i,i+2 



-2a 2 eh 2s 



0. 



2 2 

Next, observing that for all X in the space .M(j + i)2(lR) x .M(j+i)2 
II (Cw,a - I) (X) || = \\(W- \I)X + X(A- \I) || 



< 



\w-\i\ 



\A-\I\ 



\x\ 



it results that 



\£w,a-i\\ < \\w--i\ 



< C{a)h 2s . 



(13) 



Consequently, the Lyapunov endomorphism Cw,a converges uniformly to the 
identity I as h goes towards and I = o(h 2+s ) with s > 0. Using Lemma [4.11 
the operator Cw,a is invertible for h small enough. 



5 Convergence of the Discrete Method 



6 Consistency and Stability of the Discrete Method 



The consistency of the proposed method is done by evaluating the the local 
truncation error arising from the discretization of the system 



u u - Au - v xx = 0, 
v = qu xx + u 2 . 

The principal part of the first equation is 

„. , * I 2 d 4 u h 2 fd 4 u d 4 u\ r2 d 2 (Au) 



12 <9t 4 12 {dx 4 dy 4 J dt 2 
12 a x 4 ot 2 ox 2 



(14) 
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The principal part of the local error truncation due to the second part is 



It is clear that the two operators C\ v and C? uv tend toward as I and h tend 
to 0, which ensures the consistency of the method. Furthermore, the method 
is consistent with an order 2 in time and space. 

We now proceed by proving the stability of the method by applying the Lya- 
punov criterion. A linear system C(x n+ i, x n , av_i, . . .) = is stable in the 
sense of Lyapunov if for any bounded initial solution x the solution x n re- 
mains bounded for all n > 0. Here, we will precisely prove the following result. 



Lemma 2 V n : The solution (U n , V n ) is bounded independently of n whenever 
the initial solution (U°,V°) is bounded. 

We will proceed by recurrence on n. Assume firstly that ||([/ ,y°)|| < rj for 
some rj positive. Using the system ( ITT]) , we obtain 



(15) 



(16) 



C WA (U n+1 ) = C S B (U n ) + b 2 RV n - Cw^U™- 1 ) - a 2 R(F n ~ 1 + F T 
V n+l = 2c 2 RU n+1 + 2R( Cl U n + c 2 U n - v ) - V n - X + 2F™. 

where B = B — 2a 2 C\R 2 . Consequently, 

||Ay,A(^ n+1 )|| < ||£ 5iB ||.||C/"||+2|6 2 |.||^|| 

+ + 2|a 2 |(||F"- 1 || + \\F n \) 

and 

||F n+1 || < 4|c 2 |.||C/ n+1 || +4(| Cl |.||^ n || + Ical-H?/"- 1 !!) 

+||y n_1 || + ||F n_1 || + ||.F n ||. 

Next, recall that, for I = o(h s+2 ) small enough, s > 0, we have 

ai = - + 2a h 2s+2 -> ~, a 2 = -a h 2s+2 0, 

6i = 1 - 2(1 - 2a)h 2s+2 -> 1, b 2 = (1 - 2a)h 2s+2 -> 0, 
ci = (1 — 2a)h~ 2 — > oo and c 2 = a h 2s+2 — > oo, 
a 2 ci = -a(l - 2a)h 2s ->■ 0. 
As a consequence, for h small enough, 

\\C Sb \\ < 2\\B\\+2\a 2 c 1 \\\R\\ 2 < 2max(|6 1 |, 2|6 2 |) +4|a 2 Ci| < 2 + 4 = 6, (18) 



(17) 
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and the following lemma deduced from (|T3|) . 

Lemma 6.1 For h small enough, it holds for all X G A^(j+i)a(R) that 
±\\X\\ < (l-C(a)h*)\\X\\ < \\£wAX)\\ < (l + C(a)h 2s )\\X\\ < 

As a result, ([TBI) yields that 

\\\U n+1 \\ < 6\\U n \\+2\\V n \\ + ^\\U n - 1 \\+2(\\F n - 1 \\ + ||F n ||). (19) 

For n = 0, this implies that 

II^H < 12\\U°\\ + 4||V°|| +3||f/- 1 || +4(||F- 1 || + (20) 
Using the discrete initial condition 

c/° = u- 1 + ^. 

Here we identify the function ip to the matrix whom coefficients are <Pj t m — 
<f(xj,y m ). We obtain 

||^ 1 ||<r ||+/||^||. (21) 

Observing that 
it results that 

\F-^\ < |^ |W | 2 + 2^ iim )|.|C/°J + Z 2 |^ >w i 2 
and consequently, 

||F- 1 ||<||f/ || 2 + 2/||^||.||L/ ||+/ 2 ||^|| 2 . (22) 

Hence, equation (120]) yields that 

WU 1 ]] < (15 + 8/||^||)||t/°|| +4||y°|| +8||F°|| +3/||(^|| +4/ 2 ||^|| 2 . (23) 

Now, the Lyapunov criterion for stability states exactly that 

\f e > 0, 3 V > s.t; \\(U°,V°)\\ < V || (U n , V n ) \\ < e, Vn > 0. (24) 

For n = 1 and || (C/ 1 , "K 1 ) || < e, we seek an r] > for which ||(f7°, V°)|| < rj. 
Indeed, using (|23|) . this means that, it suffices to find rj such that 

8r] 2 + (19 + 8l\\(p\\)v + 3%|| + 4/ 2 ||^|| 2 - e < 0. (25) 

Choosing / small enough, we obtain a discriminant 

A = (19 + 8Z|M|) 2 - 32(3/|M| + 4Z 2 |M| 2 - e) > 
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VA - (19 + 8/IMI) n „ , . 

and a zero rji = — > 0. Consequently, it suffices to choose 

rj e]0,Th[ to obtain ([25]). Finally, flU yields that {{U 1 ]] < e. 
Now, equation (IT7|) . for n = 0, implies that 

IIV 1 !! < A(l,h,ip)\\U°f + B(l,h,cp)\\U°\\ + C(l,h,ip) + lQ\c 2 \\\V°l (26) 

where 

= 3 + 32|c 2 |, 
B(l, h,ip) = A (\d\ + 8|c 2 |(2 + lip) + lip + ^ , 



and 



C(Z, fc, y?) = 2(1 + 8|c 2 |)Z 2 |M| 2 + 4/(4|c 2 | + 1)||<^| 



h 2 

Choosing \\(U°, V°)\\ < rj, it suffices to study the inequality 

A(l, h, ip)r] 2 + (B(l, h, ip) + 16|c 2 |) ij + C(l, h,<p)-e<0. (27) 
Its discriminant satisfies 

A ~ ^ (1 + 20a + |1 - 2a\f - ^ ((4 + lQa)\\ V \\h s - e) > 



for h small enough. The zero = ^ - ( B (U i y ) + 16 N ) > . As a 

A(/,/l,^) 

consequence, for 77 g]0,t/^[ we obtain < £. Finally, for 77 g]0,77 [ with 

770 = min(?7i, r)[), we obtain IK?/ 1 , V 1 )!! < e whenever ||({7°, V°)\\ < rj. Assume 
now that the (U k , V k ) is bounded for k = 1, 2, . . . , n (by £1) whenever (U°, V°) 
is bounded by rj and let e > 0. We shall prove that it is possible to choose rj 
satisfying || (U n+1 , V n+1 ) \\ < e. Indeed, from (TT9"j) . we have 



\\U n+v \\ < 19ei + 8ef. (28) 

So, one seeks, £1 for which 8£ 2 + 19£i — e < 0. The discriminant is A = 361 + 

V361 + 32e — 19 

32e, leading to one positive zeros S\ = . Then ||t/ n+1 || < £ 

16 

whenever \\(U k , V k )\\ < e l5 fc = 1, 2, . . . , n. Next, using ffTTj) and (T281 . we have 

< (4|ci|+80|c 2 | + l)£i + (32|c 2 |+2)£ 2 . (29) 
So, it suffices as previously to choose S\ such that 

(32|c 2 | + 2) e\ + (4|ci| + 80|c 2 | + 1) £l - s < 0. 

The discriminant is A = (4|ci| + 80|c 2 | + l) 2 + 4(32|c 2 | + 2)e, leading to one 

/ ^-(41^1+801^1 + 1) 

positive zeros £1 = ; — : — : : . lhen u < e whenever 

1 1 2(32|c 2 |+2) 11 11 " 

\\(U k , V k )\\ < ei, k = 1, 2, . . . , n. Next, it holds from the recurrence hypothesis 



10 



for Bq = min(ei, e[), that there exists 7] > for which ||([/ ,y°)|| < rj implies 
that \\(U k , V k )\\ < e , for k = 1,2, ...,n, which by the next induces that 
\\(U n+1 ,V n+1 )\\ < e. 



7 Numerical implementations 



The initial data becomes 

u j, m = u ( x j, Vm, to) = u (xj, y m ) , (j, m) e I, (30) 

U lm = \ (2u (x j ,y m ) + l 2 ^u Q (xj,y m ) + {vQ) xx (xj,y m ) - 2lg(x j ,y m )) . (31) 
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